clear
rng default


all_years=1940:1:1967; 
n_years=size(all_years,2);
ntypes_w=4;
ntypes_m=4;

D_31_Gumbel=zeros(n_years,1);
D_32_Gumbel=zeros(n_years,1);
D_43_Gumbel=zeros(n_years,1);
D_21_Gumbel=zeros(n_years,1);
D_41_Gumbel=zeros(n_years,1);
D_42_Gumbel=zeros(n_years,1);

load D_R1_M
D_31_R1=D_31_M;
D_32_R1=D_32_M;
D_43_R1=D_43_M;
D_21_R1=D_21_M;
D_41_R1=D_41_M;
D_42_R1=D_42_M;

load D_R2_M
D_31_R2=D_31_M;
D_32_R2=D_32_M;
D_43_R2=D_43_M;
D_21_R2=D_21_M;
D_41_R2=D_41_M;
D_42_R2=D_42_M;



for j=1:n_years
    
year=all_years(j);

% Probability distribution of Ychosen|X
condprobsmen_temp=load('probs_allyears_data.mat',['condprobsmen' num2str(year) '_140']);
condprobsmen=condprobsmen_temp.(['condprobsmen' num2str(year) '_140']);
PYX_cond=[condprobsmen(:,2:end) condprobsmen(:,1)]; %4x5
%Note: PYX_cond(x,y)=Probability to marry woman y born in any year conditional on being man x born in 1940

% Probability distribution of Xchosen|Y
condprobswomen_temp=load('probs_allyears_data.mat',['condprobswomen' num2str(year+1) '_140']);
condprobswomen=condprobswomen_temp.(['condprobswomen' num2str(year+1) '_140']);
PXY_cond=[condprobswomen(:,2:end) condprobswomen(:,1)];  %4x5
%Note: PXY_cond(y,x)=Probability to marry man x born in any year conditional on being woman y born in 1941

% Probability distribution of X
probsmen_temp=load('probs_allyears_data.mat',['probsmen' num2str(year) '_140']);
PX=probsmen_temp.(['probsmen' num2str(year) '_140']);
%PX(x)=Probability of being man x conditional on being born in 1940

% Probability distribution of Y
probswomen_temp=load('probs_allyears_data.mat',['probswomen' num2str(year+1) '_140']);
PY=probswomen_temp.(['probswomen' num2str(year+1) '_140']);
%PY(y)=Probability of being woman y conditional on being born in 1941

[D_Gumbel,~]=param_Gumbel_closedform_MW(PYX_cond, PXY_cond, ntypes_m, ntypes_w);

D_31_Gumbel(j)=D_Gumbel(1,3);
D_32_Gumbel(j)=D_Gumbel(2,3);
D_43_Gumbel(j)=D_Gumbel(3,4);
D_21_Gumbel(j)=D_Gumbel(1,2);
D_41_Gumbel(j)=D_Gumbel(1,4);
D_42_Gumbel(j)=D_Gumbel(2,4);

end


%%
D_31_min=D_31_R1(:,1);
D_31_min(D_31_min==-Inf)=[]; %eliminate lower bounds equl to -Inf
if isempty(D_31_min) %if there are no finite lower bounds
    D_31_min=min(D_31_Gumbel)-2; %take the Gumbel estimate as y-axis lower bound
else %if there are finite lower bounds
    D_31_min=min(D_31_min)-2;  %take the lowest finite lower bound as y-axis lower bound
end
D_31_max=D_31_R1(:,2);
D_31_max(D_31_max==Inf)=[];
if isempty(D_31_max)
    D_31_max=max(D_31_Gumbel)+2;
else
    D_31_max=max(D_31_max)+2;
end

D_32_min=D_32_R1(:,1);
D_32_min(D_32_min==-Inf)=[];
if isempty(D_32_min)
    D_32_min=min(D_32_Gumbel)-2;
else
    D_32_min=min(D_32_min)-2;
end
D_32_max=D_32_R1(:,2);
D_32_max(D_32_max==Inf)=[];
if isempty(D_32_max)
    D_32_max=max(D_32_Gumbel)+2;
else
    D_32_max=max(D_32_max)+2;
end

D_43_min=D_43_R1(:,1);
D_43_min(D_43_min==-Inf)=[];
if isempty(D_43_min)
    D_43_min=min(D_43_Gumbel)-2;
else
    D_43_min=min(D_43_min)-2;
end
D_43_max=D_43_R1(:,2);
D_43_max(D_43_max==Inf)=[];
if isempty(D_43_max)
    D_43_max=max(D_43_Gumbel)+2;
else
    D_43_max=max(D_43_max)+2;
end

D_21_min=D_21_R1(:,1);
D_21_min(D_21_min==-Inf)=[];
if isempty(D_21_min)
    D_21_min=min(D_21_Gumbel)-2;
else
    D_21_min=min(D_21_min)-2;
end
D_21_max=D_21_R1(:,2);
D_21_max(D_21_max==Inf)=[];
if isempty(D_21_max)
    D_21_max=max(D_21_Gumbel)+2;
else
    D_21_max=max(D_21_max)+2;
end


D_41_min=D_41_R1(:,1);
D_41_min(D_41_min==-Inf)=[];
if isempty(D_41_min)
    D_41_min=min(D_41_Gumbel)-2;
else
    D_41_min=min(D_41_min)-2;
end
D_41_max=D_41_R1(:,2);
D_41_max(D_41_max==Inf)=[];
if isempty(D_41_max)
    D_41_max=max(D_41_Gumbel)+2;
else
    D_41_max=max(D_41_max)+2;
end

D_42_min=D_42_R1(:,1);
D_42_min(D_42_min==-Inf)=[];
if isempty(D_42_min)
    D_42_min=min(D_42_Gumbel)-2;
else
    D_42_min=min(D_42_min)-2;
end
D_42_max=D_42_R1(:,2);
D_42_max(D_42_max==Inf)=[];
if isempty(D_42_max)
    D_42_max=max(D_42_Gumbel)+2;
else
    D_42_max=max(D_42_max)+2;
end


y_axis_min=-4; %round(min([D_31_min,D_32_min, D_43_min,D_21_min, D_41_min, D_42_min]));
y_axis_max=20; %round(max([D_31_max,D_32_max, D_43_max,D_21_max, D_41_max, D_42_max]));

%% 31

for h=1:n_years
    if  D_31_R1(h,1)==-Inf
        D_31_R1(h,1)=y_axis_min;
    end
    if  D_31_R2(h,1)==-Inf
        D_31_R2(h,1)=y_axis_min;
    end
end

for h=1:n_years
    if  D_31_R1(h,2)==Inf
        D_31_R1(h,2)=y_axis_max;
    end
    if  D_31_R2(h,2)==Inf
        D_31_R2(h,2)=y_axis_max;
    end
end




figure
plot(all_years, D_31_R1(:,1), 'Color', [0 1 1], 'LineWidth',3)
hold on
upline=D_31_R1(:,2).';
polyX = [all_years, fliplr(all_years)];
polyY = [D_31_R1(:,1).', fliplr(upline)];
fill(polyX, polyY,  [0 1 1],'LineStyle','none')
hold on
%%%
upline=D_31_R2(:,2).';
polyX = [all_years, fliplr(all_years)];
polyY = [D_31_R2(:,1).', fliplr(upline)];
viu=fill(polyX,polyY,[0 0 0]);
viu.FaceAlpha = 0.2;
[dotsX,dotsY] = ndgrid(...
    min(polyX): 0.2 : max(polyX), ...
    min(polyY): 0.2 : max(polyY));      
[in,on] = inpolygon(dotsX,dotsY,polyX,polyY); 
inX = dotsX(in|on); 
inY = dotsY(in|on); 
hold on
plot(inX,inY,'k.','MarkerSize',1,'MarkerFaceColor','k'); 
hold on
%%%
plot(all_years, D_31_Gumbel, 'b', 'LineWidth',3)
hold on
%%%
yline(0,'r--','LineWidth',2)
hold off
%%%
box on
xlim([1940 1967])
ylim([ y_axis_min   y_axis_max])
set(gca,'ytick',y_axis_min:4:y_axis_max, 'FontSize',15)
xlabel('Cohort','FontSize', 15) 
ylabel(sprintf('Men contribution to D_{33,11}(\\Phi)'),'FontSize', 15)
saveas(gcf, 'D_31_M.jpg')


%% 32

for h=1:n_years
    if  D_32_R1(h,1)==-Inf
        D_32_R1(h,1)=y_axis_min;
    end
    if  D_32_R2(h,1)==-Inf
        D_32_R2(h,1)=y_axis_min;
    end
end

for h=1:n_years
    if  D_32_R1(h,2)==Inf
        D_32_R1(h,2)=y_axis_max;
    end
    if  D_32_R2(h,2)==Inf
        D_32_R2(h,2)=y_axis_max;
    end
end


figure
plot(all_years, D_32_R1(:,1), 'Color', [0 1 1], 'LineWidth',3)
hold on
upline=D_32_R1(:,2).';
polyX = [all_years, fliplr(all_years)];
polyY = [D_32_R1(:,1).', fliplr(upline)];
fill(polyX, polyY,  [0 1 1],'LineStyle','none')
hold on
%%%
upline=D_32_R2(:,2).';
polyX = [all_years, fliplr(all_years)];
polyY = [D_32_R2(:,1).', fliplr(upline)];
viu=fill(polyX,polyY,[0 0 0]);
viu.FaceAlpha = 0.2;
[dotsX,dotsY] = ndgrid(...
    min(polyX): 0.2 : max(polyX), ...
    min(polyY): 0.2 : max(polyY));      
[in,on] = inpolygon(dotsX,dotsY,polyX,polyY); 
inX = dotsX(in|on); 
inY = dotsY(in|on); 
hold on
plot(inX,inY,'k.','MarkerSize',2,'MarkerFaceColor','k'); 
hold on
%%%
plot(all_years, D_32_Gumbel, 'b', 'LineWidth',3)
hold on
%%%
yline(0,'r--','LineWidth',2)
hold off
%%%
box on
xlim([1940 1967])
ylim([ y_axis_min   y_axis_max])
set(gca,'ytick',y_axis_min:4:y_axis_max, 'FontSize',15)
xlabel('Cohort','FontSize', 15) 
ylabel(sprintf('Men contribution to D_{33,22}(\\Phi)'),'FontSize', 15)
saveas(gcf, 'D_32_M.jpg')


%% 43

for h=1:n_years
    if  D_43_R1(h,1)==-Inf
        D_43_R1(h,1)=y_axis_min;
    end
    if  D_43_R2(h,1)==-Inf
        D_43_R2(h,1)=y_axis_min;
    end
end


for h=1:n_years
    if  D_43_R1(h,2)==Inf
        D_43_R1(h,2)=y_axis_max;
    end
    if  D_43_R2(h,2)==Inf
        D_43_R2(h,2)=y_axis_max;
    end
end

figure
plot(all_years, D_43_R1(:,1), 'Color', [0 1 1], 'LineWidth',3)
hold on
upline=D_43_R1(:,2).';
polyX = [all_years, fliplr(all_years)];
polyY = [D_43_R1(:,1).', fliplr(upline)];
fill(polyX, polyY,  [0 1 1],'LineStyle','none')
hold on
%%%
upline=D_43_R2(:,2).';
polyX = [all_years, fliplr(all_years)];
polyY = [D_43_R2(:,1).', fliplr(upline)];
viu=fill(polyX,polyY,[0 0 0]);
viu.FaceAlpha = 0.2;
[dotsX,dotsY] = ndgrid(...
    min(polyX): 0.2 : max(polyX), ...
    min(polyY): 0.2 : max(polyY));      
[in,on] = inpolygon(dotsX,dotsY,polyX,polyY); 
inX = dotsX(in|on); 
inY = dotsY(in|on); 
hold on
plot(inX,inY,'k.','MarkerSize',2,'MarkerFaceColor','k'); 
hold on
%%%
plot(all_years, D_43_Gumbel, 'b', 'LineWidth',3)
hold on
%%%
yline(0,'r--','LineWidth',2)
hold off
%%%
box on
xlim([1940 1967])
ylim([ y_axis_min   y_axis_max])
set(gca,'ytick',y_axis_min:4:y_axis_max, 'FontSize',15)
xlabel('Cohort','FontSize', 15) 
ylabel(sprintf('Men contribution to D_{44,33}(\\Phi)'),'FontSize', 15)
saveas(gcf, 'D_43_M.jpg')


%% 21

for h=1:n_years
    if  D_21_R1(h,1)==-Inf
        D_21_R1(h,1)=y_axis_min;
    end
    if  D_21_R2(h,1)==-Inf
        D_21_R2(h,1)=y_axis_min;
    end
end

for h=1:n_years
    if  D_21_R1(h,2)==Inf
        D_21_R1(h,2)=y_axis_max;
    end
    if  D_21_R2(h,2)==Inf
        D_21_R2(h,2)=y_axis_max;
    end
end


figure
plot(all_years, D_21_R1(:,1), 'Color', [0 1 1], 'LineWidth',3)
hold on
upline=D_21_R1(:,2).';
polyX = [all_years, fliplr(all_years)];
polyY = [D_21_R1(:,1).', fliplr(upline)];
fill(polyX, polyY,  [0 1 1],'LineStyle','none')
hold on
%%%
upline=D_21_R2(:,2).';
polyX = [all_years, fliplr(all_years)];
polyY = [D_21_R2(:,1).', fliplr(upline)];
viu=fill(polyX,polyY,[0 0 0]);
viu.FaceAlpha = 0.2;
[dotsX,dotsY] = ndgrid(...
    min(polyX): 0.2 : max(polyX), ...
    min(polyY): 0.2 : max(polyY));      
[in,on] = inpolygon(dotsX,dotsY,polyX,polyY); 
inX = dotsX(in|on); 
inY = dotsY(in|on); 
hold on
plot(inX,inY,'k.','MarkerSize',2,'MarkerFaceColor','k'); 
hold on
%%%
plot(all_years, D_21_Gumbel, 'b', 'LineWidth',3)
hold on
%%%
yline(0,'r--','LineWidth',2)
hold off
%%%
box on
xlim([1940 1967])
ylim([ y_axis_min   y_axis_max])
set(gca,'ytick',y_axis_min:4:y_axis_max, 'FontSize',15)
xlabel('Cohort','FontSize', 15) 
ylabel(sprintf('Men contribution to D_{22,11}(\\Phi)'),'FontSize', 15)
saveas(gcf, 'D_21_M.jpg')



%% 41

for h=1:n_years
    if  D_41_R1(h,1)==-Inf
        D_41_R1(h,1)=y_axis_min;
    end
    if  D_41_R2(h,1)==-Inf
        D_41_R2(h,1)=y_axis_min;
    end
end


for h=1:n_years
    if  D_41_R1(h,2)==Inf
        D_41_R1(h,2)=y_axis_max;
    end
    if  D_41_R2(h,2)==Inf
        D_41_R2(h,2)=y_axis_max;
    end
end


figure
plot(all_years, D_41_R1(:,1), 'Color', [0 1 1], 'LineWidth',3)
hold on
upline=D_41_R1(:,2).';
polyX = [all_years, fliplr(all_years)];
polyY = [D_41_R1(:,1).', fliplr(upline)];
fill(polyX, polyY,  [0 1 1],'LineStyle','none')
hold on
%%%
upline=D_41_R2(:,2).';
polyX = [all_years, fliplr(all_years)];
polyY = [D_41_R2(:,1).', fliplr(upline)];
viu=fill(polyX,polyY,[0 0 0]);
viu.FaceAlpha = 0.2;
[dotsX,dotsY] = ndgrid(...
    min(polyX): 0.2 : max(polyX), ...
    min(polyY): 0.2 : max(polyY));      
[in,on] = inpolygon(dotsX,dotsY,polyX,polyY); 
inX = dotsX(in|on); 
inY = dotsY(in|on); 
hold on
plot(inX,inY,'k.','MarkerSize',2,'MarkerFaceColor','k'); 
hold on
%%%
plot(all_years, D_41_Gumbel, 'b', 'LineWidth',3)
hold on
%%%
yline(0,'r--','LineWidth',2)
hold off
%%%
box on
xlim([1940 1967])
ylim([ y_axis_min   y_axis_max])
set(gca,'ytick',y_axis_min:4:y_axis_max, 'FontSize',15)
xlabel('Cohort','FontSize', 15) 
ylabel(sprintf('Men contribution to D_{44,11}(\\Phi)'),'FontSize', 15)
saveas(gcf, 'D_41_M.jpg')

%% 42

for h=1:n_years
    if  D_42_R1(h,1)==-Inf
        D_42_R1(h,1)=y_axis_min;
    end
    if  D_42_R2(h,1)==-Inf
        D_42_R2(h,1)=y_axis_min;
    end
end


for h=1:n_years
    if  D_42_R1(h,2)==Inf
        D_42_R1(h,2)=y_axis_max;
    end
    if  D_42_R2(h,2)==Inf
        D_42_R2(h,2)=y_axis_max;
    end
end


figure
plot(all_years, D_42_R1(:,1), 'Color', [0 1 1], 'LineWidth',3)
hold on
upline= D_42_R1(:,2).';
polyX = [all_years, fliplr(all_years)];
polyY = [D_42_R1(:,1).', fliplr(upline)];
fill(polyX, polyY,  [0 1 1],'LineStyle','none')
hold on
%%%
upline=D_42_R2(:,2).';
polyX = [all_years, fliplr(all_years)];
polyY = [D_42_R2(:,1).', fliplr(upline)];
viu=fill(polyX,polyY,[0 0 0]);
viu.FaceAlpha = 0.2;
[dotsX,dotsY] = ndgrid(...
    min(polyX): 0.2 : max(polyX), ...
    min(polyY): 0.2 : max(polyY));      
[in,on] = inpolygon(dotsX,dotsY,polyX,polyY); 
inX = dotsX(in|on); 
inY = dotsY(in|on); 
hold on
plot(inX,inY,'k.','MarkerSize',2,'MarkerFaceColor','k'); 
hold on
%%%
plot(all_years, D_42_Gumbel, 'b', 'LineWidth',3)
hold on
%%%
yline(0,'r--','LineWidth',2)
hold off
%%%
box on
xlim([1940 1967])
ylim([ y_axis_min   y_axis_max])
set(gca,'ytick',y_axis_min:4:y_axis_max, 'FontSize',15)
xlabel('Cohort','FontSize', 15) 
ylabel(sprintf('Men contribution to D_{44,22}(\\Phi)'),'FontSize', 15)
saveas(gcf, 'D_42_M.jpg')

close all